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Abstract. Based on numerical simulations of quasistatic deformation of model granular materials, two rheological regimes 
are distinguished, according to whether macroscopic strains merely reflect microscopic material strains within the grains in 
their contact regions (type I strains), or result from instabilities and contact network rearrangements at the microscopic level 
(type II strains). We discuss the occurrence of regimes I and II in simulations of model materials made of disks (2D) or 
spheres (3D). The transition from regime I to regime II in monotonic tests such as triaxial compression is different from both 
the elastic limit and from the yield threshold. The distinction between both types of response is shown to be crucial for the 
sensitivity to contact-level mechanics, the relevant variables and scales to be considered in micromechanical approaches, the 
energy balance and the possible occurrence of macroscopic instabilities 
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INTRODUCTION 
The quasistatic limit, the rigid limit and the macroscopic limit 

Although they are modeled, at the macroscopic level, with constitutive laws in which physical time and inertia play 
no part [1,2], granular materials are most often investigated at the microscopic level by "discrete element" numerical 
methods (DEM) in which the motion of the solid bodies is determined through integration of dynamical equations 
involving masses and accelerations. Fully quasistatic approaches, in which the system evolution in configuration 
space, as some loading parameter is varied, is regarded as a continuous set of mechanical equilibrium states, are 
quite rare in the numerical literature [3, 4, 5]. It is regarded as a natural starting point, on the other hand, to perform 
suitable averages of the mechanical response of the elements of a contact network to derive the macroscopic material 
response [6]. Whether and in which cases it is possible to dispense with dynamical ingredients of the model at the 
granular level and how the quasistatic limit is approached are fundamental issues that still need clarification. 

Another set of open questions are related to the role of particle deformability. Most DEM studies include contact 
elasticity in the numerical model. Experimentally, elastic behavior is routinely measured in quasistatic tests [7] and 
sound propagation. Yet, most often, contact deflections are quite negligible in comparison with grain diameters. In 
the "contact dynamics" method [8, 9], which is used to simulate quasistatic granular rheology [10, 11, 12], grains are 
modeled as rigid, undeformable solid bodies. The influence of contact deformability on the macroscopic behavior, the 
existence of a well-defined rigid limit are thus other basic issues calling for further investigations. 

Small granular samples, as the ones used in DEM studies, often exhibit quite noisy mechanical properties. The 
approach to a macroscopic behavior expressed with smooth stress-strain curve might seem problematic, especially in 
the presence of rearrangement events, associated with instabilities at the microscopic level [13, 14]. 



The origins of strain 

The present communication shows how one may shed light on the interplay between the quasistatic, rigid and 
macroscopic limits on distinguishing two different rheological regimes and delineating their conditions of occurrence, 
in simple model materials. Macroscopic strain in solidlike granular materials has two obvious physical origins: first, 
grains deform near their contacts, where stresses concentrate (so that one models the grain interaction with a point 
force); then, grain packs rearrange as contact networks break, and then repair in a different stable configuration. We 
refer here respectively to the two different kinds of strains as type I and II. The present paper, based on numerical 



simulations of simple materials, identifies the regimes, denoted as I and II accordingly, within which one mechanism 
or the other dominates, and discusses the consequences on the quasistatic rheology of granular materials. 



NUMERICAL MODEL MATERIALS AND SIMULATION PROCEDURES 

Two sets of numerical simulation results are exploited below. Two-dimensional (2D) assemblies of polydisperse disks, 
as in Refs. [15, 3, 16], are subjected to fully stress-controlled biaxial tests, for which a quasistatic computation 
method [15, 3, 17] is exploited, in addition to standard DEM simulations. The behavior of three-dimensional (3D) 
packs of monosized spherical particles, as in Refs. [18, 19, 20] is studied in simulated triaxial compression tests, with 
special attention to strains in the quasistatic limit. Part of the results are presented in the references (mostly in some 
conference proceedings) cited just above, pending the publication of a more comprehensive study. 



Two-dimensional material and stress-controlled tests 

2D systems are simulated in order to investigate basic rheophysical mechanisms with good accuracy, in the simplest 
conceivable, yet representative, model material. Samples made of polydisperse disks in 2D, with a uniform diameter 
distribution between 0.5a and a, are first assembled on isotropically compressing frictionless particles, thus producing 
dense packs with solid fraction 4> = 0.8434 ± 3 x 10~ 4 and coordination number z close to 4 in the large system 
limit. Those values are extrapolated from data averaged on sets of samples with N = 1024, 3025 and 4900 disks. The 
samples are enclosed in a rectangular cell framed by solid walls, 2 of which are mobile orthogonally to their direction, 
which enables us to carry out biaxial compression tests (Fig. 1). Finite system effects on 4> and z are mainly due to the 
surrounding walls and can be eliminated (they are proportional to perimeter to area ratio). 




FIGURE 1. Schematic representation of the biaxial tests simulated on 2D disk samples. 02 ; 
initial isotropic pressure P, while tfi = F\ jL\ is stepwise increased. 



= F2/L2 is kept constant, equal to the 



Stress-increment controlled DEM simulations 



Once prepared in mechanical equilibrium under an isotropic pressure P, disk samples, in which contacts are now 
regarded as frictional, with friction coefficient jti = 0.25, are subjected to biaxial tests as sketched in Fig. 1. Strains 
£1 = —AL2/L2, £2 = — AZ4/L1, "volumetric" strain £ v = £] + £2 are measured in equilibrium configurations, while 
the stress deviator q = 02 — <5\ is the control parameter. We use soil mechanics conventions, for which compressive 
stresses and shrinking strains are positive, q is stepwise increased by small intervals Aq — 10~ 3 P. The contact 
model is the standard (Cundall-Strack [21]) one with normal (Kn) and tangential (Kt) stiffness constants such 
that K N = 2K T = 10 5 P. A normal viscous force is also introduced in the contact, in order to reach equilibrium 



configurations faster. After each deviator step, one waits for the next equilibrium configuration, in which forces and 
moments are balanced with good accuracy (with a tolerance below 10~ 5 Pa for forces on grains, below 10~ 5 PL for 
forces on walls). We refer to this procedure as stress-increment controlled (SIC) DEM. 



FIGURE 2. Normal (left) and tangential (right) contact behavior in 2D disk samples, as schematized with rheological elements: 
springs with stiffness constants Ajv, Kj, dashpot with damping constant T)n, plastic slider with threshold related to normal force by 
coefficient )X. 

The static elastoplastic method (hereafter referred to as SEM), amounts to dealing with the initial sample configu- 
ration as a network of springs and plastic sliders corresponding to contact behavior, as in Fig. 2 - with the dashpots 
ignored, as they play no role in statics. The evolution of the system under varying load is determined as a continuous 
trajectory in configuration space, each point of which is an equilibrium state. It has been implemented in [15, 3], and 
a similar approach was used in [4]. The algorithm will not be described here, as it is presented in [17]. It relies on 
resolution of linear system of equations, with the form of the matrix (the elastoplastic stiffness matrix) depending on 
contact status (nonsliding, sliding, open). The bases of the approach are also discussed in [5]. 

SEM calculations are possible as long as only type I strains are obtained, and the results reported here [15] 
correspond to the deviator interval < q < q\ in biaxial compressions from the chosen initial state (in which all 
tangential forces are equal to zero), in which a type I response is obtained. 



Triaxial compression tests of assemblies of N — 4000 single-sized spherical beads of diameter a are simulated 
by DEM, with the more standard procedure in which the axial strain rate e a is kept constant (hereafter strain-rate 
controlled or SRC DEM). The deviator stress, q, is measured, as a function of axial strain e a = £\, as q = <j\ — 03, 
where 0\ is the major ("axial") principal stress conjugate to e a , while the other two (lateral) principal stresses 02 = 03 
are kept equal to the initial isotropic pressure P. To allow for comparisons with laboratory experiments, the beads are 
attributed the elastic properties of glass (Young modulus E = 70 GPa, Poisson ratio v = 0.3) and friction coefficient 
/I = 0.3. The contact law is a somewhat simplified version of the Hertz-Mindlin ones [22], as in Ref. [19], which 
might be consulted for more details. It leads to favorable comparisons of elastic moduli [20] obtained in simulations 
and measured in experiments on glass beads. Preparation of cuboidal samples with periodic boundaries in all three 
directions (and thus statistically homogeneous and devoid of wall effects) under prescribed pressures in the range 
10 kPa < P < 1 MPa is detailed in Refs. [19, 23]. It is shown [19] that one may obtain, depending on the assembling 
procedure, for densities close to the random close packing limit 4> ~ 0.64 under low pressure, coordination numbers 
ranging from z — 4 (or z* — 4.5 if the "rattlers", i.e., grains carrying no force, are excluded from the count) to z = 6 
(more exactly z* — 6, with 1 or 2% of rattlers) in the limit of P — > 0. As in [19], the low-coordination systems (z* ~ 4.5) 
are referred to as "C samples" in the sequel, while those with z* — 6 are called "A samples". We can thus assess the 
influence of initial coordination number on the small strain (pre-peak) behavior of a dense material. 



The stricly quasistatic approach: SEM calculations 




Triaxial compression of 3D bead assemblies 



Dimensionless control parameters 



The contact law and th e simu lated mechanical te st lea d to the definition of useful dimensionless numbers. The 
inertia parameter I = £ a \JmjaP (in 3D) or / = t a \frnjP (in 2D) characterizes the importance of inertial effects in 
strain-rate controlled tests underpressure P (m is the grain mass). The parameter / has been used repeatedly to describe 
the state of granular materials in steady flow, both in experiments [24] and in simulations [25, 26, 27], or the departure 
from equilibrium in a slow compression [23]. It also plays a central role in the recent formulation of a successful 
constitutive law for dense granular flows [28] . 

The importance of contact deflections, relative to grain diameter a, is expressed by the stiffness number, k, which 

/ E \ 2/3 

is defined as K = Kn/P in 2D models with linear contact elasticity, and as K = ^— 2yp J W ^ ^ ^ leac ' s an( ^ 

Hertzian contacts. In both cases, typical contact deflections h satisfy h/a °< K ^ [19]. 

The three limits mentioned in the introduction can be defined as / — > (quasistatic limit), K — > °° (rigid limit), 
N — > °o (macroscopic limit). 



SIMULATION RESULTS 



Biaxial tests in 2D 



Type I response interval, quasistatic approach 
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FIGURE 3. q (normalized by P) versus e a in SIC tests on 2 samples of 3025 disks, showing a very stiff increase (confused with 
vertical axis), and then a staircase regime. Inset: detail of very small strains, with comparison of SIC DEM and SEM calculations. 

In Fig. 3, q(e a ) curves as obtained by SIC DEM are shown for two samples of 3025 disks. The curves first exhibit a 
very sharp increase of deviator q, which, as revealed once the strain scale is blown up by a factor of 10 4 in the insert, 
is in fact an interval of type I response: direct SEM calculation are possible, and coincide with DEM results. The 
smoothness of the stress variations versus strain in that range is characteristic of a continuous trajectory of equilibrium 
states in configuration space. Beyond the transition to type II strain regimes, a staircase-shaped deviator curve (Fig. 3) 
is observed, exhibiting intervals of stability (nearly vertical parts of curve in Fig. 3), separated by rearrangement events 
(horizontal parts of curve in Fig. 3) in which the system gains kinetic energy before a new stable contact network is 
formed. We could check that the SEM procedure is able to reproduce the stability intervals obtained with SIC DEM. 
On reversing the load (stepwise decreasing q), a considerably larger deviator range is accessible to SEM calculations, 
and thus in regime I, as illustrated by the two (quasi-vertical) dotted lines on the main plot of Fig. 3. 




TABLE 1. Average and standard deviation of q\ as ob- 
tained over 26 samples with N=1024, 10 samples with 
N=3025 and 6 samples with N=4900. 



| SEM 


N=l024 


N=3025 


N=490() > 


| <9l> | 


0.750 ±0.050 


0.774± 0.033 


0.786± 0.024 | 



Role of contact stiffness 

As the system, in regime I, is equivalent to a network of springs and plastic sliders (Fig. 2), type I strains are all 
inversely proportional to stiffness level k, provided the compression that decreases K does not significantly affect the 
sample geometry. The curves pertaining to different K values coincide if expressed with stress ratios and variables K£, 
as shown in Fig 4. 
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FIGURE 4. Stress ratio q/p (left) and rescaled volumetric strain ke v (right) vs. rescaled axial strain K£\ for K = 10 5 and K = 10 4 . 



Approach to the macroscopic limit 

The staircase-shaped loading curves in regime II should approach in the large sample limit a smooth stress-strain 
curve, as observed in very slow laboratory tests. To check for the approach of such a macroscopic limit, the average 
(q(£i)) and the standard deviation a(q)(ei) are computed as functions of axial strain for sets of samples of three 
different sizes, and the region of the £\ - q plane corresponding to {q){£\) — o(q) < q < (q)(£i) + o{q) is shaded 
on Fig. 5, darker zones corresponding to larger N. Fluctuations about the average curve decrease as the system size 
increases, and the insert shows that the standard deviation, as averaged over the interval < £\ < 0.02, regresses 
as TY -1 / 2 . Thus staircase curves get smoothed in the large system limit, which implies that the "stairs" become 
increasingly small and numerous: as N increases rearrangement events (microscopic instabilities) become smaller 
and smaller, but more and more frequent. A similar regression is observed for the volumetric strain curve. 

Unlike the small type I response intervals observed within regime II, the stability range q < q\ of the initial, isotropic 
structure does not dwindle as the system size increases. As shown in Table 1, the initial regime I deviator interval even 
increases a little, approaching a finite limit as N — > °°. 

Our implementation of SEM involves no creation of new contacts (although this could be taken into account in a 
more refined version). This approximation becomes exact in the rigid limit of K — > °°, because a finite strain increment 
is necessary for additional contacts to close, while type I strains scale as jc _1 . The near coincidence of SEM and DEM 
approaches, the latter involving contact creations, shows that new contacts are indeed negligible for k = 10 s . q\ thus 
represents the maximum deviator stress supported by the initial contact network, beyond which [5], due to contact 
sliding and opening, an instability or a "floppy mode" appears. The hallmark of such instabilities is the negativity of 
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FIGURE 5. Main plot: sample to sample average of q versus E\. Shaded regions extend to one standard deviation about the 
average, with, in this order, darker and darker shades of gray for N=1024, 3025, 4900. Insert: regression of fluctuations, proportional 
to N~ l l 2 . Average standard deviations a{q) and (j(e v ) over interval < £i < 0.02. 

the second-order work [5, 17], viz. 

A 2 W(AU) =AU-g-AU 

for some direction of displacement increment vector AU. Vector AU comprises all increments of grain displacements 
and rotations, and K is the stiffness matrix, which depends, via the status of contacts, on the direction of AU. 
A 2 W (AU) < implies that the increment of contact forces resulting from a small perturbation AU will accelerate 
the resulting motion, whence a spontaneous increase of kinetic energy. 

Transition stress q\ and the " critical yield analysis" approach 

One may wonder whether q\ marks the upper bound q u of the deviator interval for which contact forces balancing 
the external load (i. e., statically admissible) and satisfying Coulomb's inequality (i. e., plastically admissible) can 
be found in the network. This is the "critical yield analysis" approach to failure in structural mechanics. It is known 
that q\ and q u would coincide if the sliding in contacts where friction is fully mobilized implied dilatancy, with an 
angle equal to the friction angle (the discrete analog of an "associated" flow rule). Fig. 6 shows that q\ is well below 
q u . With a dilatant sliding rule, the material response in biaxial compression would be stiffer, and initially (rather 
paradoxically) more contractant, and the deviator would reach 1 3P (instead of about 0.8P) before failure of the initial 
contact network. 



Evolution of microscopic state variables 

As mentioned above, contact creation is negligible in regime I and the fabric evolution is essentially due to contacts 
opening, mostly in the direction of extension (direction 2). As the initial coordination number is maximal, because of 
the absence of friction in the assembling process, very few contacts are gained in the direction of compression. In the 
initial state, all contacts only bear normal force components. Friction mobilization is gradual, but the proportion of 
sliding contacts, as shown in Fig. 7 steadily increases from zero in regime I, and reaches an apparent plateau in regime 
II. This means that the interval of elastic response is, strictly speaking, reduced to naught, even though the stress- 
strain curve can approximately be described as elastic in a very small range. The appearance of sliding contacts can 
effectively reduce the degee of static indeterminacy in the system. If the status is assumed to be fixed for all contacts, 
the Coulomb condition, satisfied as an equality, reduces the number of independent contact components from 2dN c (in 
2D) to 2N C — Xs, with % s the number of sliding contacts among a total of N c . It has been speculated [29, 30] that failing 
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FIGURE 6. Comparison of SEM calculation with the normal (parallel, curves marked "n. a." for non associated) and the 
"associated" (dilatant, curves marked "a.") sliding rule in contacts in sample with 1024 disks. 
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FIGURE 7. In a sample with N=3025, degree of force indeterminacy h (solid line) and proportion of sliding intergranular contacts 
Xs/N c (dotted line), versus axial strain E\. 

contact networks (regime II) should correspond to vanishing force indeterminacy. The data of Fig. 7 provide evidence 
against such a prediction, as h stabilizes to about 600, a moderate (10% of the total number of degrees of freedom in a 
sample of 3025 disks), yet finite value. 

3D triaxial tests 

The simulations reported here compare dense states A (high coordination number) and C (low coordination number). 
State A is similar to the dense disk sample studied in the previous section, as both were initially assembled with 
frictionless grains. (A samples, once packed under low pressure, are nevertheless compressed to the desired confining 
pressure with the value fi = 0.3 of the friction coefficient used in the triaxial tests [19]). Pressure values correspond 
to glass beads, and vary between 10 kPa (k = 39000) and 1 MPa (k = 1800). We first check for the approach of 
the quasistatic and the macroscopic limit in 3D, strain-rate controlled DEM simulations, then discuss the influence of 
coordination number, and regimes I and II, in the light of the previous 2D study. 





Reproducibility, quasistatic limit 



As the system size increases, sample to sample fluctuations should regress, as checked in 2D (Fig. 5). Our 3D results 
are based on 5 samples of 4000 beads of each type, and Fig. 8 checks for stress-strain curve reproducibility in both A 
and C cases, for small axial strains. Thanks to the fully periodic boundary conditions [19], the macroscopic mechanical 
behavior is quite well defined with N = 4000. The approach to the quasistatic limit, in SRC tests can be assessed on 
checking for the innocuousness of the dynamical parameters, i.e., inertial number /, and reduced damping parameter 
£. £ is defined as the ratio of the viscous damping constant in a contact to its critical level, given the instantaneous 
value of the stiffness constant. We found it convenient to use a constant value of £ in our simulations, as in [19]. Fig. 8 
also shows that provided inertial number /, characterizing dynamical effects, is small enough, both / and £ become 
irrelevant. Fig. 8 shows that the quasistatic limit is correctly approached for / < 10 -3 , quite a satisfactory result, given 
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FIGURE 8. Left: small strain part of q(e a ) curves for 5 different samples of each type, A (top curves) and C (bottom ones) with 
N = 4000 beads. Right: q{e a ) and E v (e a ) curves in one type C sample for the different values of £ and / indicated. 

that usual laboratory tests with t a ~ 10~ 5 s _1 correspond to / < 10~ 8 . 



Influence of initial coordination number 



Fig. 9 compares the behavior of initial states A and C, in triaxial compression with P = 100 kPa (k ~ 6000). 
Although, conforming to the traditional view that the peak deviator stress is determined by the initial sample density, 
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FIGURE 9. q(e a ) (left scale) and e v (e a ) (right scale) curves for A and C states under P = 100 kPa. Averages over 5 samples of 
4000 spherical grains. 

maximum q values are very nearly identical in systems A and C, the mobilization of internal friction is much more 



gradual for C. For A, the initial rise of deviator q for small axial strain is quite steep, and the volumetric strain variation 
becomes dilatant almost immediately, for e a ~ 10~ 3 . In [20] it was shown that measurements of elastic moduli provide 
information on coordination numbers. It is thus conceivable to infer the rate of deviator increase as a function of axial 
strain from very small strain (<~ 10~ 5 or below [7, 20]) elasticity. Most experimental curves obtained on sands, which 
do not exhibit q maxima or dilatancy before e a ~ 0.01, are closer to C ones. However, some measurements on glass 
bead samples [31] do show fast rises of q at small strains, somewhat intermediate between numerical results of types 
A and C. 



Influence of contact stiffness 

The small strain (say e a < 5.10~ 4 ) interval for A samples, with its fast q increase, is in regime I, as one might 
expect from 2D results on disks. This is readily checked on changing the confining pressure. Fig. 10 shows the 
curves for triaxial compressions at different P values (separated by a factor VTO) from 10 kPa to 1 MPa, with a 
rescaling of the strains by the stiffness parameter K, in one A sample. Their coincidence for q/P < 1 evidences a 
wide deviator range in regime I. For larger strains, curves separate on this scale, and tend to collapse together if 
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FIGURE 10. Left: q(e a )/P and E v (e a ) curves for one A sample and different P values. Strains on scale {P/Po) 2 ^ °= k: -1 , 
Pq — 100 kPa. Right: q(e a )/P for the same P values in one C sample. Inset: detail with blown-up e scale, straight lines corresponding 
to Young moduli in isotropic state. 

q/P, e v are simply plotted versus e a . The strain dependence on stress ratio is independent from contact stiffness. This 
different sensitivity to pressure is characteristic of regime II. Fig 10 also shows that it applies to C samples almost 
throughout the investigated range, down to small deviators (a behavior closer to usual experimental results than type 
A configurations). At the origin (close to the initial isotropic state, see inset on fig. 10, right plot), the tangent to the 
curve is given by the elastic (Young) modulus of the granular material, E m , and therefore q/P scales with K, but curves 
quickly depart from this behaviour (around q = 0.2P). The approximately elastic range [20] is quite small, as observed 
in experiments [7, 32, 33]. 



Calculations with a fixed contact list 

Within regime I, the mechanical properties of the material can be successfully predicted on studying the response of 
one given set of contacts. Those might slide or open, but the very few new contacts that are created can be neglected. 
To check this in simulations, one may restrict at each time step the search for interacting grains to the list of initially 
contacting pairs. Fig. 1 1 compares such a procedure to the complete calculation. The curve marked "NCC" for no 
contact creation is indistinguishable from the other one for q > 0.8. We thus check that, in regime I, the macroscopic 
behavior is essentially determined by the response of a fixed contact network. 
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FIGURE 11. Very small strain part of q(£ a ) curve in one A sample, showing beginning of unloading curves (arrows). Curve 
marked NCC was obtained on calculating the evolution of the same sample without any contact creation. 



Fig. 1 1 also shows that the small strain response of A samples, within regime I, close to the initial state, is already 
irreversible: type I strains are not elastic. An approximately elastic behavior is only observed for very small strains, 
as depicted in the inset of Fig. 10 (right part). In this small interval near the initial equilibrium configuration, the 
stress-strain curve is close to its initial tangent, defined by the elastic modulus. Moduli [20] can be calculated from the 
stiffness matrix of contact networks. One may also check that the unloading curves shown on Fig. 11 (and the ones 
of Fig. 3 in 2D as well) comprise a small, approximately elastic part, with the relevant elastic modulus (the Young 
modulus for a triaxial test at constant lateral stress) defining the initial slope. At the microscopic level, a small elastic 
response is retrieved upon reversing the loading direction because contacts stop sliding. The elastic range is strictly 
included in the larger range of type I behavior. 



Finally, let us note that regimes I and II also differ by the importance of sample to sample fluctuations: curves in 
Fig. 8 (left plot) pertaining to the different samples of type A or C are confused as long as q < I. IP (case A) or 
q < 0.3P (case C), which roughly corresponds to the transition from regime I to regime II. Larger fluctuations imply 
that the characteristic length scale associated with the displacement field (correlation length) is larger in regime II. 
Whether and in what sense rearrangements triggered by instabilities in regime II, in a material close to the rigid limit 
(large fc), can be regarded as local events is still an open issue. 



Numerical studies thus reveal that the two regimes, in which the origins of strain differ, exhibit contrasting properties. 
Although the reported studies in 2D and 3D differ in many respects (linear versus Hertzian contacts, wall versus 
periodic boundaries, SIC versus SRC DEM), the same phenomena were observed in both cases. Regime I corresponds 
to the stability range of a given contact structure. It is larger in highly coordinated systems. It is observed in the 
beginning of monotonic loading tests, in which the deviator stress increases from an initial isotropic configuration, and 
also after changes in the direction of load increments (hence a loss in friction mobilization). Strains, for a given stress 
level, are then inversely proportional to contact stiffnesses. The deviator range in regime I, q < q\, in usual monotonic 
tests, is stricly larger than the small elastic range, but strictly smaller than the maximum deviator. It does not vanish 
in the limit of large systems, unlike in the singular case of rigid, frictionless particle assemblies [13, 34]. Regime I is 
limited by the occurrence of elastoplastic instabilities in the contact network and does not coincide with the prediction 
of the critical yield approach. In regime I, the work of the externally applied load is constantly balanced by the one of 
contact forces, so that the kinetic energy approaches zero in the limit of slow loading rates. A remarkable consequence 



Type I strains and elastic reponse 



Fluctuations and length scale 



CONCLUSION 



is that the instability condition based on the negativity of the macroscopic second-order work [35] is never fulfilled, 
as macroscopic and microscopic works coincide, and the latter is positive. In regime II, network rearrangements are 
triggered by instabilities and some bursts of kinetic energy are observed [ 14] . Larger fluctuations witness longer-ranged 
correlations in the displacements. The microscopic origin of macroscopic strains, which are independent on contact 
elasticity for usual stiffness levels K, lies in the geometry of grain packings. 

On attempting to predict a macroscopic mechanical response from packing geometry and contact laws, the infor- 
mation about which kind of strain should dominate is crucial. 

A promising perspective is the study of correlated motions associated with rearrangement events. 
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